---
title: ANOVAs for a two-condition ERP experiment 
output: 
  html_notebook
fontsize: 11pt
geometry: margin=1in
---
This code assumes that you have already exported data from EEGLAB/ERPLAB. I've been using the ERP measurement tool in ERPLAB to extract peak amplitudes and mean amplitudes. This code assumes you have an output file from the measurement tool, and that it is in long format.

#Preprocessing

Let's load some libraries. 
Tidyverse loads the Wickham libraries dplyr, tidyr, ggplot2.
EZ gives us ANOVA functions that have the Greenhous-Geisser sphericity correction, which we will need for factors with more than two levels.
Knitr gives us lots of things, including the kable() function for displaying data frames in a nice looking table.

```{r}
library(tidyverse)
library(ez)
library(knitr)
```

Now let's read in the data. This assumes that the data is in your working directory. We also bind all of the data sets together into one.

I use ERPLAB's ERP measurement tool to create a dataset in long format that contains the mean amplitude between 300ms and 500ms. It would also be possible to use R to calculate these things... but since ERPLAB can do it, it seems easier for now to just use ERPLAB to generate these. This script is wri

```{r}
d = read.table("300-500.txt", header=TRUE)

```

The next step is to edit the dataset with factors that will allow us to run the ANOVAs that we want to run. This means choosing the channels that we want to analyze, creating subsets like midline vs medial vs lateral, marking hemispheres (left/right), and electrode positions. The chunks below do this.

First we label channels yes/no based on whether we want to analyze them. This arrangement is close to the Neville et al. 1991 set of electrodes. This is one of the steps where the researcher exercises researcher degrees of freedom to reduce the dimensionality of the data.

```{r}
Nchannels = rep(c("no",
                  "yes",
                  "yes",
                  "yes",
                  "no",
                  "no",
                  "no",
                  "yes",
                  "yes",
                  "no",
                  "no",
                  "no",
                  "yes",
                  "yes",
                  "yes",
                  "no",
                  "no",
                  "no",
                  "yes",
                  "yes",
                  "no",
                  "no",
                  "no",
                  "yes",
                  "yes",
                  "yes",
                  "no",
                  "no",
                  "no",
                  "yes",
                  "yes",
                  "no"), times=length(levels(factor((d$bini)))) * length(levels(factor(d$ERPset))))
```

Next we divide the electrodes into subsets. One common way to do this is to divide into three groups: midline (Xz), medial (X3/4), and lateral (X7/8). There are other ways. This is one of the steps where the researcher exercises researcher degrees of freedom to reduce the dimensionality of the data.

```{r}
locationAnalysis = rep(c("no",
                         "midline",
                         "medial",
                         "lateral",
                         "no",
                         "no",
                         "no",
                         "medial",
                         "lateral",
                         "no",
                         "no",
                         "no",
                         "midline",
                         "medial",
                         "lateral",
                         "no",
                         "no",
                         "no",
                         "medial",
                         "lateral",
                         "no",
                         "no",
                         "no",
                         "midline",
                         "medial",
                         "lateral",
                         "no",
                         "no",
                         "no",
                         "medial",
                         "lateral",
                         "no"), times=length(levels(factor((d$bini)))) * length(levels(factor(d$ERPset))))
```

Next we can divide the non-midline channels by hemisphere (by creating a hemisphere factor, we can test for laterality in the effect)

```{r}
hemisphere = rep(c("no",
                   "no",
                   "left",
                   "left",
                   "no",
                   "no",
                   "no",
                   "left",
                   "left",
                   "no",
                   "no",
                   "no",
                   "no",
                   "left",
                   "left",
                   "no",
                   "no",
                   "no",
                   "right",
                   "right",
                   "no",
                   "no",
                   "no",
                   "no",
                   "right",
                   "right",
                   "no",
                   "no",
                   "no",
                   "right",
                   "right",
                   "no"), times=length(levels(factor((d$bini)))) * length(levels(factor(d$ERPset))))
```

Finally, we can divide electrodes by their anterior/posterior position. This lets us test for anteriority. I just use numbers (1 is anterior, higher is posterior)

```{r}
electrodePosition = rep(c(
  "no",
  "1",
  "1",
  "1",
  "no",
  "no",
  "no",
  "2",
  "2",
  "no",
  "no",
  "no",
  "3",
  "3",
  "3",
  "no",
  "no",
  "no",
  "3",
  "3",
  "no",
  "no",
  "no",
  "2",
  "2",
  "2",
  "no",
  "no",
  "no",
  "1",
  "1",
  "no"), times=length(levels(factor((d$bini)))) * length(levels(factor(d$ERPset))))
```

Now we just add these vectors to the dataframe, and do some minor editing to the column names.

```{r}

d.full = data.frame(d, Nchannels, locationAnalysis, hemisphere, electrodePosition) 

#change the name of columns ----
colnames(d.full)[4]="condition"
colnames(d.full)[5]="subject"
colnames(d.full)[1]="voltage"

#convert to factors if not already ----
d.full$condition = factor(d.full$condition)
d.full$subject = factor(d.full$subject)
```

The first step before running the ANOVAs is to choose the subset of channels. There is only one set of options here, but one could imagine having different sets of channels.

```{r}
#subset by channels ----
d.N = d.full %>%
    filter(Nchannels=="yes") %>%
    droplevels
```


#Option 1: Run a separate ANOVA on each subset of electrodes (midline, medial, lateral)

This is what Neville et al. 1991 do. With three ANOVAs, we probably should use Dunn correction or some other form of correction, but this is unfortunately rare in the literature. 

###midline electrodes

```{r}
mid = d.N %>%
  filter(locationAnalysis=="midline") %>%
  droplevels

mid.omni = mid %>%
  filter(condition == "1" | condition=="2") %>%
  droplevels

anova.mid = ezANOVA(mid.omni, dv=voltage, wid=subject, within=.(condition, electrodePosition))

```

###medial electrodes

```{r}
medial = d.N %>%
  filter(locationAnalysis=="medial") %>%
  droplevels()

medial.omni = medial %>%
  filter(condition == "1" | condition=="2") %>%
  droplevels()

anova.medial = ezANOVA(medial.omni, dv=voltage, wid=subject, within=.(condition, hemisphere, electrodePosition))

```

###lateral electrodes

```{r}
lateral = d.N %>%
  filter(locationAnalysis=="lateral") %>%
  droplevels()

lateral.omni = lateral %>%
  filter(condition == "1" | condition=="2") %>%
  droplevels()

anova.lateral = ezANOVA(lateral.omni, dv=voltage, wid=subject, within=.(condition, hemisphere, electrodePosition))

```


#Option 2: Collapse electrodes together into groups, and run an ANOVA with a factor called electrode that has levels for each group.

Here we will make 6 clusters of electrodes, hemisphere (left, right) x electrodePosition (1, 2, 3 - anterior, central, posterior). We probably could do this through some dplyr magic, but I am just going to hardcode it by creating subsets.

d1 = filter(d.N, hemisphere == "left" & electrodePosition == "1")
d2 = filter(d.N, hemisphere == "left" & electrodePosition == "2")
d3 = filter(d.N, hemisphere == "left" & electrodePosition == "3")
d4 = filter(d.N, hemisphere == "right" & electrodePosition == "1")
d5 = filter(d.N, hemisphere == "right" & electrodePosition == "2")
d6 = filter(d.N, hemisphere == "right" & electrodePosition == "3")

left.anterior = d1 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
left.anterior$electrode = rep("left anterior", times=nrow(left.anterior))

left.central = d2 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
left.central$electrode = rep("left central", times=nrow(left.central))

left.posterior = d3 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
left.posterior$electrode = rep("left posterior", times=nrow(left.posterior))

right.anterior = d4 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
right.anterior$electrode = rep("right anterior", times=nrow(right.anterior))

right.central = d5 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
right.central$electrode = rep("right central", times=nrow(right.central))

right.posterior = d6 %>%
  group_by(condition, subject) %>%
  summarize(voltage = mean(voltage))
right.posterior$electrode = rep("right posterior", times=nrow(right.posterior))

##append the 6 together
d.electrodes = rbind(left.anterior, left.central, left.posterior, right.anterior, right.central, right.posterior)
d.electrodes$electrode = factor(d.electrodes$electrode)

###midline electrodes

```{r}

anova.electrode = ezANOVA(d.electrodes, dv=voltage, wid=subject, within=.(condition, electrode))

```
